%==============================================================================
% This code computes equilibrium with CBDC under different CBDC rate
% Cournot competition in both deposit and loan markets
% with a perfectly competitive interbank market and
% banks free entry.
%  COPY RIGHT: Yu Zhu
%==============================================================================


clear;close all
options = optimset('maxiter',5e5,'tolX',1e-10,'tolfun',1e-10,'maxfunevals',5e5);

%==========================================
%Define parameters
%==========================================
name = 'LN87_08_final';
load(['pars_post_' name '.mat']);
name = [name,'_free_entry'];
com_opt = 1; % loan market competition 1 if cournot 0 if perfect competition
%year_option = 0; % 1: using banking data from 2007 and 2008. 0: using banking data from 2014-2019
location='figures/'

alpha1  =   pars(4);
alpha2  =   pars(5);
alpha3  =   pars(6);
beta    =   pars(7);
epsilon =   0.001;
sigma   =   pars(1);
cost    =   pars(8);                   
B       =   pars(2);
Nb      =    pars(3); %number of banks
A       =   pars(10);  %productivity
eta     =   pars(9);                                       %bank's cost to manage deposits
theta   =   pars(11);
rq_ratio=   pars(12); %reserve ratio
i_reserve = pars(13);
year_option = pars(14)
%% Define the name of the saved files: last two number indicates reserve requirement ratio and the interest on reserves
if year_option==1
    pi_m = 0.03346;pi_h=pi_m;
    name = [name,'_', num2str(floor(rq_ratio*100)),'_',num2str(floor(i_reserve*100)),'_', num2str(beta*100), '_07_08'];
else
    pi_m = 0.01515;pi_h=pi_m;
    name = [name,'_', num2str(floor(rq_ratio*100)),'_',num2str(floor(i_reserve*100)),'_', num2str(beta*100)];
end
nrhogrid=  10000;
i_m     = 1/beta*(1+pi_m)-1;
i_h     = 1/beta*(1+pi_h)-1;
pi_eff  = (1+pi_m)/(1+i_reserve)-1;




%==============================================
%Define functions
%===============================================

u = @(x) ((x+epsilon).^(1-sigma)-epsilon.^(1-sigma))./(1-sigma); %DM utility
du = @(x)(x+epsilon).^(-sigma); %derivative of DM utility
d2u=@(x)(x+epsilon).^(-sigma-1)*(-sigma);
c=@(x)x;                                                         % DM cost
dc=@(x)1;%   derivative of DM cost
d2c=@(x)0;

y_star = fminsearch(@(x)(du(x)-dc(x)).^2,1);                     % y_star
g_func=@(x) (1-theta).*u(x)+theta.*c(x);                         %trading proctocol 
dg_func=@(x) (1-theta).*du(x)+theta.*dc(x);
d2g_func=@(x)(1-theta).*d2u(x)+theta.*d2c(x);
y_grid = 0:0.001:y_star;
p_grid = g_func(y_grid);
p_star = p_grid(end);
y_func = @(p)interp1(p_grid,y_grid,min(p,max(p_grid)));          %trading quantity 
p_func=@(x)min(x,max(p_grid));                                   %trading payment
lambda_func = @(z) max(du(y_func(z))./dg_func(y_func(z))-1,0);          % liquidity premium

f=@(k)A*k.^(eta);%production
if com_opt==1
    df_inv = @(rho,Nb)(A*eta*(1-(1-eta)/Nb)./(1+rho)).^(1/(1-eta));
    df=@(k,Nb)A*eta*(1-(1-eta)/Nb)*k.^(eta-1); %marginal production
else
    df_inv = @(rho,Nb)(A*eta./(1+rho)).^(1/(1-eta));
    df=@(k,Nb)A*eta*k.^(eta-1); %marginal production
end

dlambda_func=@(z)((d2u(y_func(z)).*dg_func(y_func(z))-du(y_func(z)).*d2g_func(y_func(z)))./dg_func(y_func(z)).^3).*(z<=p_grid(end));
dinv=@(z)beta*(alpha2*lambda_func(z)+1);
ddinv=@(z)beta*alpha2*dlambda_func(z);
z_grid = linspace(0,p_grid(end),100000);
lambda_grid = lambda_func(z_grid);
gamma_grid=linspace(0,20,nrhogrid);
x0=[0,0,0];
x0_R=x0;
opt=1;
%icbdc_grid=linspace(-0.01,(1+pi_m)/beta-1,1000);
icbdc_grid=linspace(0,0.02,1000);

picbdc=pi_m;
D_grid = linspace(0,max(p_grid),5e3);

pass_variables

ygrid = D_inv_demand(i_m,D_grid,func_var);
rho_grid = linspace(-.05,rho_bar,100);
[LS,LD,rho_grid] = LSLD(rho_grid,func_var,pi_m,Nb,ygrid);
LD_cournot =( A*eta*((eta-1)/Nb+1)./(1+rho_grid)).^(1/(1-eta));
% eq_plot=figure
% hold on
% h1=plot(1+rho_grid,LS,'black','linewidth',2)
% plot([1+rho_grid(end),1+rho_grid(end)],[LS(end),LS(end)*2],'black','linewidth',2)
% h2=plot(1+rho_grid, LD,'blue','linewidth',2)
% %plot(1+rho_grid, LD_cournot,'linewidth',2)
% 
% ylim([0,1.7])
% xlabel('$R_{I}$','interpreter','latex')
% ylabel('$L$','interpreter','latex')
% legend([h1,h2],'Loan Supply','Loan Demand','location','northwest')
% figure_halfpage
% saveas(eq_plot,[location,'eq_plot_entry'],'epsc')
% saveas(eq_plot,[location,'eq_plot_entry'],'png')

%====================================================================
%Compute entry cost
%====================================================================
 y=SS_eq_noCBDC(func_var,pi_eff,Nb,ygrid);

psi= beta*(alpha2*lambda_func(y(2))+alpha3*lambda_func(y(1)+y(2)))+beta;
drate=1/psi-1;
loan= fzero(@(x)df(x,Nb)-1-y(end),[1e-10,100]);

eq_final_nocbdc1=[y,drate,loan];

 y=SS_eq_noCBDC(func_var,pi_eff,Nb+1,ygrid);

psi= beta*(alpha2*lambda_func(y(2))+alpha3*lambda_func(y(1)+y(2)))+beta;
drate=1/psi-1;
loan= fzero(@(x)df(x,Nb+1)-1-y(end),[1e-10,100]);

eq_final_nocbdc2=[y,drate,loan];
rl_rate1 = (1+eq_final_nocbdc1(:,3))/(1-(1-eta)/Nb)-1;
rl_rate2 = (1+eq_final_nocbdc2(:,3))/(1-(1-eta)/(Nb+1))-1;
welfare_banker1 = max(1+rl_rate1-1/(1+pi_eff),0).*eq_final_nocbdc1(:,end)-eq_final_nocbdc1(:,2)...
    -(cost-(1/(1+pi_eff))).*eq_final_nocbdc1(:,2)./(1+eq_final_nocbdc1(:,4));

welfare_banker2 = max(1+rl_rate2-1/(1+pi_eff),0).*eq_final_nocbdc2(:,end)-eq_final_nocbdc2(:,2)...
    -(cost-(1/(1+pi_eff))).*eq_final_nocbdc2(:,2)./(1+eq_final_nocbdc2(:,4));

welfare_banker1=welfare_banker1/Nb;
welfare_banker2=welfare_banker2/(Nb+1);
entry_cost = (welfare_banker1+welfare_banker2)/2; %entry cost defined as the middle point of the profits with Nb and Nb+1 banks.
entry_cost =welfare_banker1;

for i=1:length(icbdc_grid)
    i
    
    if i==1
    y=SS_eq_noCBDC(func_var,pi_eff,Nb,ygrid);
    end
    eq_picbdc=(1+picbdc)/(1+icbdc_grid(i))-1;
    y_reserve = SS_eq_CBDCR(func_var,pi_m, i_reserve, icbdc_grid(i),Nb,ygrid);
    
    i_m=(1+pi_m)/beta-1;
    psi= beta*(alpha2*lambda_func(y(2))+alpha3*lambda_func(y(1)+y(2)))+beta;
    drate(i)=1/psi-1;
    loan= fzero(@(x)df(x,Nb)-1-y(end),[1e-10,100]);
    
    eq_final_nocbdc(i,:)=[y,drate(i),loan];
    x0=y;
    x0_R=y_reserve;
    psi_reserve =  beta*(alpha2*lambda_func(y_reserve(2))+alpha3*lambda_func(y_reserve(1)+y_reserve(2)))+beta;
    drate_reserve(i)=1/psi_reserve-1;
    loan_reserve=fzero(@(x)df(x,Nb)-1-y_reserve(end),[1e-10,100]);
    
    %===============================================
    %Compute equilibrium where CBDC is not a reserve
    %===============================================
    icbdc_real  = (1+picbdc)/(icbdc_grid(i)+1)/beta-1;
    if icbdc_real>psi/beta-1
        eq_cbdc(i,:)=[y,drate(i),loan];
        drate_cbdc(i)=drate(i);
        eliquidity(i,1) = y(2);
        rl_rate1 = (1+eq_cbdc(i,3))/(1-(1-eta)/Nb)-1;
        profit_cbdc_nb = max(1+rl_rate1-1/(1+pi_eff),0).*eq_cbdc(i,end)-eq_cbdc(i,2)...
            -(cost-(1/(1+pi_eff))).*eq_cbdc(i,2)./(1+eq_cbdc(i,4));
        profit_cbdc(i,:) = [profit_cbdc_nb,Nb];
    else
        tic
        z_eq= SS_eq( alpha1,alpha2,alpha3,z_grid,lambda_grid,i_m,icbdc_real);
        nb = Nb;
        lsupply_max=(1-rq_ratio)*(1+icbdc_real)*beta*z_eq(1,2);
        rho_hat=max(1/(1+pi_eff)-1, ((icbdc_grid(i)+1)/(1+picbdc)-rq_ratio/(1+pi_eff)+cost)/(1-rq_ratio)-1);
            ld_low= fzero(@(x)df(x,nb)-1-rho_hat,[1e-10,100]);
            drate_cbdc(i)=(icbdc_grid(i)+1)/(1+picbdc)-1;
            toc
            if ld_low<lsupply_max
                deposit =  ld_low/(icbdc_real+1)/beta/(1-rq_ratio);
                eq_cbdc(i,:)=[z_eq(1),deposit,rho_hat,drate_cbdc(i),ld_low];
            else
                rhocbdc = df(lsupply_max,nb)-1;
                
                eq_cbdc(i,:)=[z_eq,rhocbdc,drate_cbdc(i),lsupply_max];
                
            end
            eliquidity(i,1) = z_eq(2);
            rl_rate1 = (1+eq_cbdc(i,3))/(1-(1-eta)/nb)-1;
            profit_cbdc_nb = max(1+rl_rate1-1/(1+pi_eff),0).*eq_cbdc(i,end)-eq_cbdc(i,2)...
                -(cost-(1/(1+pi_eff))).*eq_cbdc(i,2)./(1+eq_cbdc(i,4));
            profit_cbdc(i,:) = [profit_cbdc_nb,nb]; 
    
            
            
        

    end
    
    %===============================================
    %Compute equilibrium where CBDC is a reserve
    %===============================================
    if icbdc_real>psi_reserve/beta-1
        eq_cbdc_R(i,:)=[y_reserve,drate_reserve(i),loan_reserve];
        drate_cbdc_R(i)=drate_reserve(i);
        eliquidityR(i,1) = y_reserve(2);
        rl_rate1 = (1+eq_cbdc_R(i,3))/(1-(1-eta)/Nb)-1;
         profit_cbdc_nbR = max(1+rl_rate1-max(1/(1+pi_eff),...
                (1+icbdc_grid(i))/(1+pi_m)),0).*eq_cbdc_R(i,end)-eq_cbdc_R(i,2)...
                -(cost-max(1/(1+pi_eff),(1+icbdc_grid(i))/(1+pi_m))).*eq_cbdc_R(i,2)./(1+eq_cbdc_R(i,4));
        profit_cbdcR(i,:) = [profit_cbdc_nbR,Nb];
    else
        tic
        z_eq= SS_eq( alpha1,alpha2,alpha3,z_grid,lambda_grid,i_m,icbdc_real);
        
        lsupply_maxR=(1-rq_ratio)*(1+icbdc_real)*beta*z_eq(1,2);
        rho_hatR=max(1/(1+pi_eff)-1, ((icbdc_grid(i)+1)/(1+picbdc)-rq_ratio*max(1/(1+pi_eff),1/(1+eq_picbdc))+cost)/(1-rq_ratio)-1);
         nb = Nb
            
            ld_lowR= fzero(@(x)df(x,nb)-1-rho_hatR,[1e-10,100]);
            drate_cbdc_R(i)=(icbdc_grid(i)+1)/(1+picbdc)-1;
            toc
            if ld_lowR<lsupply_maxR
                depositR =  ld_lowR/(icbdc_real+1)/beta/(1-rq_ratio);
                eq_cbdc_R(i,:)=[z_eq(1),depositR,rho_hatR,drate_cbdc_R(i),ld_lowR];
                
                %break
            else
                rhocbdcR = df(lsupply_maxR,nb)-1;
                
                eq_cbdc_R(i,:)=[z_eq,rhocbdcR,drate_cbdc_R(i),lsupply_maxR];
                
            end
            eliquidityR(i,1) = z_eq(2);
            rl_rate1 = (1+eq_cbdc_R(i,3))/(1-(1-eta)/nb)-1;
            profit_cbdc_nbR = max(1+rl_rate1-max(1/(1+pi_eff),...
                (1+icbdc_grid(i))/(1+pi_m)),0).*eq_cbdc_R(i,end)-eq_cbdc_R(i,2)...
                -(cost-max(1/(1+pi_eff),(1+icbdc_grid(i))/(1+pi_m))).*eq_cbdc_R(i,2)./(1+eq_cbdc_R(i,4));
                profit_cbdcR(i,:) = [profit_cbdc_nbR,nb];

        

    end
    
end
markup=(alpha1*min(eq_cbdc(1,1),max(p_grid))/y_func(eq_cbdc(1,1)) +alpha2*min(eq_cbdc(1,2),max(p_grid))/y_func(eq_cbdc(1,2))...
    +alpha3*min(eq_cbdc(1,2)+eq_cbdc(1,1),max(p_grid))/y_func(eq_cbdc(1,2)+eq_cbdc(1,1)))/(alpha1+alpha2+alpha3);
    
Output_CBDC = alpha1*eq_cbdc(:,1)+alpha2*eliquidity+alpha3*min(eq_cbdc(:,1)+eliquidity,p_star)+B+eq_cbdc(:,end)+f(eq_cbdc(:,end))-eq_cbdc(:,2)...
    +(eq_cbdc(:,2)./(1+eq_cbdc(:,4))-eq_cbdc(:,end))*(1+i_reserve)/(1+pi_m);
Output_CBDCR = alpha1*eq_cbdc_R(:,1)+alpha2*eliquidityR+alpha3*min(eq_cbdc_R(:,1)+eliquidityR,p_star)+B+eq_cbdc_R(:,end)+f(eq_cbdc_R(:,end))-eq_cbdc_R(:,2)...
     +(eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4))-eq_cbdc_R(:,end)).*(1+max(i_reserve,icbdc_grid(:)))/(1+pi_m);


output_index=(Output_CBDC/Output_CBDC(1)-1)*100;
output_indexR=(Output_CBDCR/Output_CBDCR(1)-1)*100;
macro_implication=figure
psi_D = eq_cbdc(:,2)./(1+eq_cbdc(:,4));
psi_D_R = eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4));
subplot(3,3,1)
hold on
h1=plot(icbdc_grid*100,(psi_D./psi_D(1)-1)*100,'linewidth',2);
h2=plot(icbdc_grid*100,(psi_D_R./psi_D_R(1)-1)*100,'--red','linewidth',2);
xlim([-1,3])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Deposit')

subplot(3,3,2)
hold on
plot(icbdc_grid*100,(eq_cbdc(:,end)./eq_cbdc(1,end)-1)*100,'linewidth',2)
plot(icbdc_grid*100,(eq_cbdc_R(:,end)./eq_cbdc(1,end)-1)*100,'--red','linewidth',2)
xlim([-1,3])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Loan')

subplot(3,3,3)
hold on
plot(icbdc_grid*100,output_index,'linewidth',2)
plot(icbdc_grid*100,output_indexR,'--red','linewidth',2)
xlim([-1,3])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Output')

ndrate=((1+eq_cbdc(:,4))*(1+pi_m)-1)*100;
ndrateR=((1+eq_cbdc_R(:,4))*(1+pi_m)-1)*100;
subplot(3,3,4)
hold on
plot(icbdc_grid*100,((1+eq_cbdc(:,4))*(1+pi_m)-1)*100,'linewidth',2)
plot(icbdc_grid*100,((1+eq_cbdc_R(:,4))*(1+pi_m)-1)*100,'--red','linewidth',2)
xlim([-1,3])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Deposit Rate')
if com_opt==1
    rl_rate = (1+eq_cbdc(:,3))./(1-(1-eta)./profit_cbdc(:,2))-1;
    rl_rateR = (1+eq_cbdc_R(:,3))./(1-(1-eta)./profit_cbdcR(:,2))-1;
else
    rl_rate = eq_cbdc(:,3);
    rl_rateR = eq_cbdc_R(:,3);
end

nlrate=((1+rl_rate)*(1+pi_m)-1)*100;
nlrateR=((1+rl_rateR)*(1+pi_m)-1)*100;
subplot(3,3,5)
hold on
plot(icbdc_grid*100,nlrate,'linewidth',2)
plot(icbdc_grid*100,nlrateR,'--red','linewidth',2)
xlim([-1,3])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Loan Rate')

subplot(3,3,6)
hold on
plot(icbdc_grid*100,nlrate-ndrate,'linewidth',2)
plot(icbdc_grid*100,nlrateR-ndrateR,'--red','linewidth',2)
xlim([-1,3])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Spread')

hL = subplot(3,3,8);
poshL = get(hL,'position');     % Getting its position

lgd = legend(hL,[h1,h2],'CBDC not as Reserve','CBDC as Reserve', 'NumColumns',2);
set(lgd,'position',poshL);      % Adjusting legend's position
axis(hL,'off');                 % Turning its axis off

saveas(macro_implication,[location,'interest_bearing_',name],'epsc')
saveas(macro_implication,[location,'interest_bearing_',name],'png')

%% Not as reserves
macro_implication=figure
psi_D = eq_cbdc(:,2)./(1+eq_cbdc(:,4));
psi_D_R = eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4));
subplot(2,3,4)
hold on
h1=plot(icbdc_grid*100,(psi_D./psi_D(1)-1)*100,'linewidth',2);
xlim([0,.8])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Deposit')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

subplot(2,3,5)
hold on
plot(icbdc_grid*100,(eq_cbdc(:,end)./eq_cbdc(1,end)-1)*100,'linewidth',2)
xlim([0,.8])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Loan')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

subplot(2,3,6)
hold on
plot(icbdc_grid*100,output_index,'linewidth',2)
xlim([0,.8])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Output')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

ndrate=((1+eq_cbdc(:,4))*(1+pi_m)-1)*100;
ndrateR=((1+eq_cbdc_R(:,4))*(1+pi_m)-1)*100;
subplot(2,3,1)
hold on
plot(icbdc_grid*100,((1+eq_cbdc(:,4))*(1+pi_m)-1)*100,'linewidth',2)
xlim([0,.8])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Deposit Rate')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

%nlrate=((1+eq_cbdc(:,3))*(1+pi_m)-1)*100;
%nlrateR=((1+eq_cbdc_R(:,3))*(1+pi_m)-1)*100;
subplot(2,3,2)
hold on
plot(icbdc_grid*100,nlrate,'linewidth',2)
xlim([0,.8])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Loan Rate')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

subplot(2,3,3)
hold on
plot(icbdc_grid*100,nlrate-ndrate,'linewidth',2)
xlim([0,.8])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Spread')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

saveas(macro_implication,[location,'interest_bearing_nonreserve_',name],'epsc')
saveas(macro_implication,[location,'interest_bearing_non_reserve_',name],'png')
% Buyer welfare

%alternative welfare calculation
welfare=alpha1*(u(y_func(eq_cbdc(:,1)))-p_func((eq_cbdc(:,1))))+alpha2*(u(y_func(eliquidity))-p_func((eliquidity)))...
+alpha3*(u(y_func(eq_cbdc(:,1)+eliquidity))-p_func((eq_cbdc(:,1)+eliquidity)))-(B/2-eq_cbdc(:,2)+eq_cbdc(:,2)./(1+eq_cbdc(:,4)))+B/2*log(B/2)...
-((1+i_reserve)/(1+pi_m)-1).*(eq_cbdc(:,2)./(1+eq_cbdc(:,4))-eq_cbdc(:,end))-cost*(eliquidity-eq_cbdc(:,2))./(1+eq_cbdc(:,4));

welfareR=alpha1*(u(y_func(eq_cbdc_R(:,1)))-p_func((eq_cbdc_R(:,1))))+alpha2*(u(y_func(eliquidityR))-p_func((eliquidityR)))...
+alpha3*(u(y_func(eq_cbdc_R(:,1)+eliquidityR))-p_func((eq_cbdc_R(:,1)+eliquidityR)))+B/2*log(B/2)-(B/2-eq_cbdc_R(:,2)+eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4)))...
-((1+max(i_reserve,icbdc_grid(:)))./(1+pi_m)-1).*(eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4))-eq_cbdc_R(:,end))-cost*(eliquidityR-eq_cbdc_R(:,2))./(1+eq_cbdc_R(:,4));

welfare_func=@(x)alpha1*(u(x*y_func(eq_cbdc(1,1)))-p_func((eq_cbdc(1,1))))+alpha2*(u(x*y_func(eq_cbdc(1,2)))-p_func((eq_cbdc(1,2))))...
+alpha3*(u(x*y_func(eq_cbdc(1,1)+eq_cbdc(1,2)))-p_func((eq_cbdc(1,1)+eq_cbdc(1,2))))+B/2*log(x*B/2)-(B/2-eq_cbdc(1,2)+eq_cbdc(1,2)/(1+eq_cbdc(1,4)))...
-((1+i_reserve)/(1+pi_m)-1).*(eq_cbdc(1,2)./(1+eq_cbdc(1,4))-eq_cbdc(1,end));


for i=1:length(icbdc_grid)
    welfare_comp(i) = fzero(@(x)welfare_func(x)-welfare(i),[0.5,1.5]);
    welfare_compR(i) = fzero(@(x)welfare_func(x)-welfareR(i),[0.5,1.5]);
end
h3=figure;
hold on
plot(icbdc_grid*100,(welfare_comp-1)*100,'linewidth',2)
plot(icbdc_grid*100,(welfare_compR-1)*100,'--red','linewidth',2)
xlabel('$i_e$','interpreter','latex')
xlim([0,5])
legend('CBCD not Reserve','CBDC as Reserve','location','northwest')
figure_halfpage
saveas(h3,[location,'welfare_b',name],'png')

%entrepreneur welfare
welfare_e = f(eq_cbdc(:,end))-(1+eq_cbdc(:,3)).*eq_cbdc(:,end);

welfare_eR = f(eq_cbdc_R(:,end))-(1+eq_cbdc_R(:,3)).*eq_cbdc_R(:,end);

h2=figure
hold on
plot(icbdc_grid*100,(welfare_e/welfare_e(1)-1)*100,'linewidth',2)
plot(icbdc_grid*100,(welfare_eR/welfare_eR(1)-1)*100,'--red','linewidth',2)
xlabel('$i_e$','interpreter','latex')
xlim([0,5])
[m,n]=max(welfare_e/welfare_e(1))
[m,nR]=max(welfare_eR/welfare_eR(1))
figure_halfpage
legend('CBDC not Reserve','CBDC as Reserve','location','southeast')

saveas(h2,[location,'welfare_e',name],'png')


%banker welfare
welfare_banker = max(1+rl_rate-1/(1+pi_eff),0).*eq_cbdc(:,end)-eq_cbdc(:,2)-(cost-(1/(1+pi_eff))).*eq_cbdc(:,2)./(1+eq_cbdc(:,4));
welfare_bankerR = max(1+rl_rateR-max(1/(1+pi_eff),(1+icbdc_grid(:))/(1+pi_m)),0).*eq_cbdc_R(:,end)-eq_cbdc_R(:,2)-(cost-max(1/(1+pi_eff),(1+icbdc_grid(:))/(1+pi_m))).*eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4));
h1=figure
hold on
plot(icbdc_grid*100,(welfare_banker/welfare_banker(1)-1)*100,'linewidth',2)
plot(icbdc_grid*100,(welfare_bankerR/welfare_bankerR(1)-1)*100,'--red','linewidth',2)
xlabel('$i_e$','interpreter','latex')
xlim([0,2])
figure_halfpage
legend('CBDC not Reserve','CBDC as Reserve')
saveas(h1,[location,'welfare_bankers',name],'png')

%welfare of sellers

welfareS=alpha1*(-c(y_func(eq_cbdc(:,1)))+p_func((eq_cbdc(:,1))))+alpha2*(-c(y_func(eliquidity))+p_func((eliquidity)))...
+alpha3*(-c(y_func(eq_cbdc(:,1)+eliquidity))+p_func((eq_cbdc(:,1)+eliquidity)))-B/2+B/2*log(B/2);
welfareSR=alpha1*(-c(y_func(eq_cbdc_R(:,1)))+p_func((eq_cbdc_R(:,1))))+alpha2*(-c(y_func(eliquidityR))+p_func((eliquidityR)))...
+alpha3*(-c(y_func(eq_cbdc_R(:,1)+eliquidityR))+p_func((eq_cbdc_R(:,1)+eliquidityR)))+B/2*log(B/2)-B/2;
welfare_funcS=@(x)alpha1*(-c(y_func(eq_cbdc(1,1)))+p_func((eq_cbdc(1,1))))+alpha2*(-c(y_func(eq_cbdc(1,2)))+p_func((eq_cbdc(1,2))))...
+alpha3*(-c(y_func(eq_cbdc(1,1)+eq_cbdc(1,2)))+p_func((eq_cbdc(1,1)+eq_cbdc(1,2))))+B/2*log(x*B/2)-B/2;
for i=1:length(icbdc_grid)
    i
    welfare_compS(i) = fzero(@(x)welfare_funcS(x)-welfareS(i),[0.5,1.5]);
    welfare_compSR(i) = fzero(@(x)welfare_funcS(x)-welfareSR(i),[0.5,1.5]);
end
h4=figure
hold on
plot(icbdc_grid*100,(welfare_compS-1)*100,'linewidth',2)
plot(icbdc_grid*100,(welfare_compSR-1)*100,'--red','linewidth',2)
xlabel('$i_e$','interpreter','latex')
xlim([0,5])
legend('CBCD not Reserve','CBDC as Reserve','location','northwest')
figure_halfpage
saveas(h4,[location,'welfare_S',name],'png')



%=============================
%Compute certain numbers
%=============================
[loan_max,nmax] = max(eq_cbdc(:,end)./eq_cbdc(1,end));
[loan_maxR,nmaxR] = max(eq_cbdc_R(:,end)./eq_cbdc(1,end));
i_max=icbdc_grid(nmax);
i_maxR=icbdc_grid(nmaxR);
id = find(eq_cbdc(:,end)./eq_cbdc(1,end)>=1,1,'last');
id1 = find(eq_cbdc(:,end)./eq_cbdc(1,end)>1+1e-8,1,'first');
idR1 = find(eq_cbdc_R(:,end)./eq_cbdc(1,end)>1,1,'first');

idR = find(eq_cbdc_R(:,end)./eq_cbdc(1,end)>=1,1,'last');
i_max1=icbdc_grid(id);
i_min1=icbdc_grid(id1);
i_minR1=icbdc_grid(idR1);
i_maxR1=icbdc_grid(idR);

[output_max,nmax] = max(output_index);
[outputR_max,nmaxR] = max(output_indexR);
ioutput_max=icbdc_grid(nmax);
ioutput_maxR=icbdc_grid(nmaxR);
id = find(output_index>=0,1,'last');
id1 = find(output_index>0,1,'first');
idR = find(output_indexR>=0,1,'last');
idR1 = find(output_indexR>0,1,'first');

iy_max1=icbdc_grid(id);
iy_min1=icbdc_grid(id1);
iy_maxR1=icbdc_grid(idR);
iy_minR1=icbdc_grid(idR1);


disp('=====================================================================================');
disp(['If CBDC does not serve as reserves, it increases Bank Lending if its rate is between 0.3049% and ' num2str(i_max1*100) '%']);
disp(['It increases output if its rate is between 0.3049% and ' num2str(iy_max1*100) '%']);
disp(['The CBDC rate maximizing lending is ' num2str(i_max*100) '% and that maximizing output is ' num2str(ioutput_max*100) '%' ]);
disp(['The maximum increase in lending is ' num2str((loan_max-1)*100) '% and that maximizing output is ' num2str(output_max) '%' ]);
disp('')
disp('=====================================================================================');
disp(['If CBDC serves as reserves, it increases Bank Lending if its rate is between 0.3049% and ' num2str(i_maxR1*100) '%']);
disp(['It increases output if its rate is between 0.3049% and ' num2str(iy_maxR1*100) '%']);
disp(['The CBDC rate maximizing lending is ' num2str(i_maxR*100) '% and that maximizing output is ' num2str(ioutput_maxR*100) '%' ]);
disp(['The maximum increase in lending is ' num2str((loan_maxR-1)*100) '% and that maximizing output is ' num2str(outputR_max) '%' ]);

